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Inertial particles are often observed to be trapped, temporarily or permanently, by recirculation cells which 
are ubiquitous in natural or industrial flows. In the limit of small particle inertia, determining the conditions 
of trapping is a challenging task, as it requires a large number of numerical simulations or experiments to test 
various particle sizes or densities. Here, we investigate this phenomenon analytically and numerically in the 
case of heavy particles (e.g. aerosols) at low Reynolds number, to derive a trapping criterion that can be used 
both in analytical and numerical velocity fields. The resulting criterion allows to predict the characteristics of 
trapped particles as soon as single-phase simulations of the flow are performed. Our analysis is valid for two- 
dimensional particle-laden flows in the vertical plane, in the limit where the particle inertia, the free-fall terminal 
velocity, and the flow unsteadiness can be treated as perturbations. The weak unsteadiness of the flow generally 
induces a chaotic tangle near heteroclinic or homoclinic cycles if any, leading to the apparent diffusion of fluid 
elements through the boundary of the cell. The critical particle Stokes number St c below which aerosols also 
enter and exit the cell in a complex manner has been derived analytically, in terms of the flow characteristics. It 
involves the non-dimensional curvature-weighted integral of the squared velocity of the steady fluid flow along 
the dividing streamline of the recirculation cell. When the flow is unsteady and St > St c , a regular motion 
takes place due to gravity and centrifugal effects, like in the steady case. Particles driven towards the interior 
of the cell are trapped permanently. In contrast, when the flow is unsteady and St < St c , particles wander in a 
chaotic manner in the vicinity of the border of the cell, and can escape the cell. 
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I. INTRODUCTION 

The settling of inertial particles in plane two-dimensional 
flows characterized by quasi-steady recirculation cells is a sit¬ 
uation of wide interest in natural or industrial flows. These 
cells can be due for example to the presence of obstacles of 
any kind, to temperature gradients, density gradients, or wind 
on water flows (see for example Stommel 1 , Chen & Fung 2 , 
Fung 3 , Maxey & Corrsin 4 , Cerisier et al. 5 , to cite but a few). 
They are characterized by well-defined scales and shapes, and 
their lifetime is much larger than their turnover time. They 
have been shown to play a key role in the dynamics of inertial 
particles, that is tiny objects which do not follow exactly the 
fluid motion, like aerosols, sediments or even biological ob¬ 
jects like plankton 1 . Particles which penetrate into these cells 
often perform a few rotations, then exit or get deposited on the 
walls if any. They can enter under the effect of gravity, inertia, 
or under the effect of any slight perturbation occurring while 
the particle passes near the boundary of the recirculation cell. 
This trapping process can be permanent, for example if parti¬ 
cles get deposited on walls, or temporary if particles exit the 
cell. 

In most situations, theoretical issues related to the com¬ 
plex dynamics of inertial particles can be investigated with 
a reasonable accuracy by means of numerical computations 
or analytical models (see for example the review paper by 
Cartwright et al. 1 ). The former generally requires heavy com¬ 
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puting resources, and consists in solving the Navier-Stokes 
equation (for the fluid motion) coupled to Lagrangian or Eu- 
lerian models for the dispersed phase. In most cases, such 
computations are time-consuming, even if particles do not af¬ 
fect the flow. The prediction of the degree of contamination 
of recirculation cells by inertial particles, which requires nu¬ 
merous computations with varying parameters, is therefore a 
challenging task. In this paper, we combine analytical and 
numerical methods to derive an analytical expression for the 
critical Stokes number (and therefore the critical size or mass) 
of particles which are likely to contaminate a given cell. This 
method is valid if the flow field can be written as a steady 
mean flow plus any small time-dependent perturbations. 

The work is presented in three steps. First, we consider 
the transport equations and derive an analytical criterion al¬ 
lowing to determine the conditions under which inertial par¬ 
ticles cross the separatrix, that is the limiting streamline be¬ 
tween the cell and the rest of the flow domain (section II). It 
is then used to predict particle capture in a simple analytical 
flow (section III). Finally, the criterion is applied to accurate 
numerical solutions of the Navier-Stokes equations, coupled 
to a Lagrangian particle tracking algorithm, in section IV. 


II. ASYMPTOTIC ANALYSIS OF SEPARATRIX 
CROSSING 

The structure of steady plane incompressible flows is 
the classical structure of 2D volume preserving dynam¬ 
ical systems, where the hamiltonian corresponds to the 
streamfunction . These flows are therefore composed of sad¬ 
dle or elliptic stagnation points, the former being joined by 
streamlines, and the latter being surrounded by them. The 
recirculation cells are the elliptic stagnation points and the 
closed streamlines around them. The ’’boundaries” of the cells 
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FIG. 1: Sketch of the stagnation points A and B of the steady flow 
u () , joined by a separatrix streamline (panel (i)). Vector n is the 
unit vector perpendicular to uo, such that (uo, nj is anti-clockwise 
about the z axis. Some examples which will be treated in this paper: 
circular cell (ii), cavity flow (iii), backward facing-step (iv). 


are formed by limiting streamlines joining stagnation points, 
and forming homoclinic or heteroclinic cycles. These cycles 
are separatrices playing a key role in the momentary or per¬ 
manent trapping of inertial particles. 

Let L be the typical length scale of the flow (e.g. the overall 
scale of a recirculation cell), and let U be the typical velocity 
of the flow. The non-dimensional motion equation of non- 
Brownian spherical particles with a small Reynolds number 
and a density p p much larger than the fluid density py reads 


dv 

dt 



- v) + fV 4 


(1) 


where v is the velocity of the center-of-mass of the particle, 
St = D 2 p p U/(18pL) is the Stokes number (where p, is the 
dynamic viscosity of the fluid), Fr = U 2 /{gL) is the Froude 
number, and g is the unit vector in the direction of gravity. 
The generalization of this theory to particles with a density of 
the order of that of the fluid, like sediments in water, is among 
the perspectives of this work. Also, particle/wall interactions, 
either by direct contact, electrostatic forces or hydrodynamic 
interactions, will be neglected. We will assume that particles 
can slip on the walls, if any. Throughout the paper we will 
consider that particles have a finite but small inertia, that is 


0 < St < 1. 


The flow is composed of a steady component u 0 plus a small 
unsteady component eui (with 0 < e < 1). This unsteady 
term can be due to some noise, or to ambient turbulence with 
a weak turbulence intensity. We exploit the fact that St <C 1 
and look for a solution of Eq. (1) of the form 

v = u 0 + St vi + 0(St 2 ). (2) 

In addition, we assume that even though unsteadiness, inertia 
and sedimentation are small perturbations, none of these ef¬ 
fects dominates the other two. Therefore, St is of the order of 


the amplitude of the unsteady component, that is St = 0(e), 
and the free-fall velocity St/Fr is also of the order of the 
unsteady component eui. In particular, the Froude number 
Fr is held fixed as e and St tend to zero. Approximation (2) 
manifests the fact that if particles are sufficiently small, their 
velocity v(t) is always close to the steady fluid velocity Uo(x) 
at the particle position x(f), and gravity, inertia as well as flow 
unsteadiness can be treated as a perturbations 9 . Under these 
conditions we obtain a reduced dynamical equation for iner¬ 
tial particles (see also Refs. 10,11 ): 

St 

x(t) = v(f) ~ u 0 (x) +eui(x,i) + — g- St (u 0 .V)u 0 , 

(3) 

plus quadratic combinations of e and St. This equation will 
be used below. 

We now consider a separatrix streamline (also denoted as 
dividing streamline) of the steady flow Uo joining two sad¬ 
dle stagnation points A and B, with the flow from A to B 
(Fig. 1). In the following, the separatrix will be simply de¬ 
noted by AB. The dynamics near AB can be investigated 
by using classical methods originally developed for separatrix 
maps 12-14 . We consider a solution x(f ) of the complete pertur¬ 
bative equation (3), and determine the variations of the unper¬ 
turbed streamfunction t/>o[x(f)] (defined by Uq = Vx z)) 
along this perturbed trajectory. Also, we will consider the 
’’undisturbed streamfunction seen by the particle”, that is 
the variation of 0 O at discrete times r„ and r n+1 , where r n 
and r n +1 are two consecutive times when the trajectory x(f) 
passes closest to the stagnation points A and B respectively: 

AV’o = ^o(x(r n+ i)) - V>o(x(t„)). (4) 

Using Eq. (3) this quantity reads: 

Ai/'o = J V^o-(eui(x,f) + St (u 0 .V)u 0 ) dt. 

(5) 

The sign of Aiftg contains direct information about the be¬ 
havior of the particle near the heteroclinic cycle containing 
AB: a constant sign will indicate that the particle penetrates 
(or exits) the cell once and for all, whereas oscillating signs 
correspond to a chaotic motion: particles cross the separatrix 
AB in a complex manner (see Chirikov 14 , Guckenheimer & 
Holmes 15 , Del Castillo-Negrete 20 ). 

The variation of the undisturbed streamfunction ALq is cal¬ 
culated as follows. We consider a cartesian coordinates sys¬ 
tem (x, y) in the plane of the flow, with unit vector x in the 
horizontal direction, and unit vector y directed upward (that 
is, y = — g). Let t n £]r„, r n +i [ be the time when x(t ) passes 
nearest to some reference point C on separatrix AB (with C 
not equal to A nor B). By writing that x(f) ~ q(£ — t n ), 
where q(f) is an exact solution of the unperturbed dynamics 
on AB with q(0) = C, we are led to: 

Of 

A'i/’o ~ AV >0 + ^r;( AB x g) z - St (|u 0 | 2 ) c , (6) 

where AB is the displacement vector joining A and B, and: 

r+oo 

At/)£ = ez. qxui(q (t),t + t n )dt, (7) 
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is the variation of streamfunction for fluid points, and 


for all n is: 


/ -t-oo 

qxq dt (8) 

-OO 

is related to centrifugal effects and will be discussed below. 
In these calculations we have introduced the unit vector z = 
x x y, which is perpendicular to the plane of the flow. Also, 
(AB x g). z = (AB x y). z = —AB.it, and will be 
denoted as —xab in the following. 

These results have been obtained by making use of the fact 
that, by definition of q, we have q = Uq and q = (uo.V)uo. 
Also, the dynamics of q(f) is very slow near the stagnation 
points A and B, so that integrals over [r„ — f n ,r„+i — t n ] 
have been replaced by integrals over [— oo, oo] 13 - 14 . The co¬ 
efficient (|uo| 2 ) c is the effect of the curvature of separatrix 
AB on the dynamics of inertial particles. Indeed, the local 
curvature is i?” 1 = z.q x q/|q| 3 so that (|uo| 2 ) c is exactly 
zero for straight separatrices, and is strictly non-zero for sep- 
aratrices with a constant curvature. Inertial particles, experi¬ 
encing a centrifugal force along curved separatrices, tend to 
drift towards the exterior of the recirculation cell, and this ef¬ 
fect brings a constant contribution into A'tif. The coefficient 
(|uo| 2 ) c can be re-written as a curvilinear integral along the 
separatrix: 



where s is the arc-length, and we made use of the fact that 
|q|<if = |ciq| = ds. The last integral is a curvature-weighted 
integral of the undisturbed squared fluid velocity, and this is 
why we have chosen to denote it by the bracket (-) c . Note that 
this integral is not always positive, as the curvature is a signed 
quantity: when non-zero, f? _1 is positive for anti-clockwise 
streamlines, and negative otherwise. Under these conditions 
Eq. (6) leads to: 

St 

AV >0 - A-00 - — XAB - St (|u 0 | 2 ) c . (10) 

Eq. (10) shows that flow unsteadiness, gravity and curvature 
bring additive contributions into Ai/jg . A positive A0 q means 
that particles moving in the vicinity of separatrix AB will 
drift, during the time interval [r n ,r„+i], towards the left of 
the streamline (i.e. towards n, see Fig. 1 (i)). Negative contri¬ 
butions correspond to particles drifting towards its right-hand- 
side (i.e. towards —n). The flow unsteadiness brings a varying 
sign. The contribution of gravity has the sign of — xab and, 
since particles are heavier than the fluid here, it always corre¬ 
sponds to particles driven downward. The curvature-weighted 
integral of the squared velocity (|u 0 | 2 ) c is positive (negative) 
if AB forms an anti-clockwise (clockwise) curve. In both 
cases it corresponds to a centrifugal effect. 

Finally, let ea denote the peak value of the variation of 
streamfunction of perturbed fluid points trajectories along the 
separatrix: 


ea = max \Af{(t n )\ 

—oo<t n <oo 

A sufficient condition for the sign of A'tlf to remain constant 


St > St c 


ea 

\x A B/Fr + (|u 0 | 2 ) c |' 


( 11 ) 


where a only depends on the structure of the flow and mani¬ 
fests the effect of unsteadiness. The curvature-weighted inte¬ 
gral (| Uo | 2 ) c also is independent of the particle characteristics. 
It only depends on the geometry of the separatrix AB and on 
the steady velocity distribution along this streamline. 

All variables in the results above have been set non- 
dimensional by using some flow length-scale L and velocity 
U. The dimensional counterpart of Eq. (11) is: 

Tp>Tpc= \g X AB + (\Uo\ 2 )cl ^ 


where t p = St X L/U is the response time of the iner¬ 
tial particle ; Xab is the dimensional ^-coordinate of vector 
AB ; U 0 is the dimensional undisturbed fluid velocity vec¬ 
tor ; 5 = ea x UL has the dimension of a diffusion coeffi¬ 
cient (?n 2 /s). It is the flux of fluid points within the stochastic 
zone near separatrix AB, due to the flow unsteadiness (see for 
example Beigie et al. I6,17 , Solomon et al. 18 ). Note that this 
quantity is not related to inertial particles and only depends 
on the structure and on the unsteadiness of the fluid flow. 

The criterion in Eq. (12) can also be written as D > D c , 
where D c is the critical diameter of the particle: 

/ \ 1/2 

D c = ^18 ^r pc J . (13) 

This expression will be used in the next sections, and com¬ 
pared to numerical simulation results. 

Structure of invariant manifolds. To understand further 
the physical meaning of these analytical results, consider the 
structure of the invariant manifolds of particle’s dynamics 
(sketched in Fig. 2). Assuming that particles are allowed to 
slip on the wall, separatrix AB and walls form a heteroclinic 
cycle. We consider Poincare sections of the particle dynam¬ 
ics: x(f„ + kT), where k is an integer, T is the period of some 
external perturbation and t n is any arbitrary time where the 
Poincare section starts. Because the leading-order dynamics 
in (3) has hyperbolic saddle-points A and B, the Poincare map 
of the perturbed dynamics will also display hyperbolic points 
A and B in the vicinity of the undisturbed ones, the eigenval¬ 
ues of which have the same signs than the unperturbed ones 
provided the perturbation is small enough. Therefore, an un¬ 
stable invariant manifold W u emerges from A and a stable 
manifold converges towards B. The two manifolds, which 
coincided with the separatrix streamline in the steady case, no 
longer coincide in the unsteady dynamics in general. The exis¬ 
tence of intersection points between both manifolds is related 
to the existence of zeros in A0 q (t n ) (see for example Ref. 1 ’’). 
If A'tpQ (i n ) is strictly negative for all t n , then manifolds do 
not intersect each other and take the form sketched in Fig. 
2(a) and (c) [in the case of the cavity and of the backward¬ 
facing step with AB towards positive x]. Any low-.S7 particle, 
injected sufficiently close to A, will be driven towards the in¬ 
terior of the cell and never exit. This happens in the steady 
case, or in the unsteady case with St > St c . In contrast, if 
A Uq (in) has simple zeros when t n varies, then the manifolds 
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FIG. 2: Sketch of the invariant manifolds of the dynamics of par¬ 
ticles in the cavity flow (a) — (b) and in the backward facing step 
flow (c) — (d), when particles are assumed to slip on the walls. 
In the steady case, or in the unsteady case with St > St c (that is 
D > D c ), manifolds do not intersect and take the form (a) and (c) 
if gravity dominates inertial curvature effects: the motion is regular 
and particles entering the cell cannot exit. In the unsteady case with 
St < St c (that is D < Dc), manifolds take the form ( b) and (d): 
the motion is chaotic, some particles can enter and exit the cell in 
a complex manner. Case (e) corresponds to a regular motion when 
centrifugal effect dominates. Panel (/) is a sketch of a trajectory 
wandering in the stochastic zone near the separatrix. 


will have an infinite number of intersection points (Fig. 2(b) 
and (d))\ particles injected near A can have a chaotic dynam¬ 
ics near the heteroclinic loop, and enter and exit the cell in an 
unpredictable manner. This kind of behaviour happens in the 
unsteady case when St < St c . 

Note that the regular behaviour can also happen when cen¬ 
trifugal effects dominate gravity and unsteadiness. Indeed, 
if AB is clockwise ((|u 0 | 2 ) < 0) and —(|u 0 | 2 ) dominates 
both the unsteady term A i/jq and the gravity term xab/Ft, 
then At/>g(i n ) is always positive: particles will be centrifuged 
away irrespective of their Stokes numbers, provided it is small 
enough. In this regime, sketched in Fig. 2(e), any particle 
injected close enough to the heteroclinic cycle will be cen¬ 
trifuged away. 


III. APPLICATION TO AN ANALYTICAL FIELD AND 
VALIDATION OF THE NUMERICAL APPROACH 


The method, and its applicability to numerical velocity 
fields, is first tested by means of an elementary analytical ve¬ 
locity field: 


tp(x,y,t) = ijj 0 (x - e sin tot, y) 


~ i/j 0 (x,y) + v 0 (x,y)esmut + 0(e 2 ) (14) 

where the steady (nondimensional) streamfunction tfo is taken 
to be: 


(15) 


and vo(x, y) = —dipo/dx is the corresponding vertical veloc¬ 
ity. This steady flow is sketched in Fig. 1 (ii): it consists in 
a pair of half-circular cells with a unit radius. The perturbed 
streamfunction if corresponds to a periodic rigid-body trans¬ 
lation of the flow (15) in the horizontal direction. 

We first consider the upper dividing streamline, that is the 
upper half circle AB in the half plane y > 0. Gravity is 
along — y: it pushes particles within the cell, and is opposed 
to centrifugal effects. By inserting the analytical expression of 
uq(x, y) into Eq. (9), and setting x = sin <f> and y = cos (f>, the 
curvature-weighted integral is readily calculated. We obtain 
(|uo| 2 ) c = —27 t, the negative sign being due to the fact that 
the upper dividing streamline is clockwise. Moreover, for this 
simple flow, a solution q = (q x (t), q y (t)) of the unperturbed 
dynamics on AB can be readily obtained, and we get 

Qx(t) = tanh 2t and q y (t) = ^ . (16) 

By inserting these functions into (7) we obtain an exact ex¬ 
pression for A^q , and: 

=* LoshU/4) + 2St ^ - b’ (17) 

where we have used xab = 2. Therefore, in the steady case 
(e = 0), A^o (t n ) is of the sign of (it — p-) for all t n : this 
means that, if Fr < 1 / 7 r gravity drives inertial particles inside 
cell irrespective of the (small) Stokes number St. In contrast, 
if Fr > I/ 7 r centrifugal effects dominate and particles slip on 
the arc AB and slowly drift away from it: any low-St parti¬ 
cle coming from the upstream region ’’bounces” on the cell. 
This last situation corresponds to the case where the unstable 
manifold W u emerging from A is towards the exterior of the 
cell and never intersects W s [as sketched in Fig. 2)(e) for the 
step]. Note that a peculiar case emerges if Fr = 1/ir. here 
gravity and centrifugal effect balance each other and At/»o(fn) 
is an oscillating function, whatever the (small) Stokes number 
is. 

For the lower dividing streamline AB of the circular cell, 
i.e. the half-circle AB in the plane y < 0, we have (|uo| 2 ) c = 
+27T, and we get 

A-tffa) c : /4) S>nU n -2St(n + F). (18) 

Hence, the critical Froude number does not exist here, and 
both gravity and centrifugal effects push the particles out of 
the cell, downward. 

In the unsteady case (e 7 ^ 0), formulas (17) and (18) im¬ 
ply that the undisturbed streamfunction seen by the particle 
oscillates with t n if: 


St < Stf(uj) 


7T ui 2 e 

17r ± 1/Frj 4cosh(7rw/4) ’ 


(19) 


the ± sign corresponding to the upper (+) and lower (-) sepa¬ 
ratrix respectively. The curves St+ (ui) and St~ (w) have been 
plotted in Fig. 3 when Fr = 0.1 (< 1 /it) and e = 0.05: 
they have a maximum when ui ~ 2.63, and decay rapidly 
when ui increases above this value. These curves will be com¬ 
pared to numerical results below. Finally, the theory can be 
applied to the horizontal separatrix of the cellular flow (with 


tfi 0 (x,y) = y(x 2 + y 2 - 1), 
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FIG. 3: Continuous lines : critical Stokes numbers Stt (solid line) 
and St~ (dashed line), for the circular separatrices of the analyt¬ 
ical cell flow (14), obtained from the theoretical expressions (19). 
Symbols : the critical Stokes numbers obtained from numerical sim¬ 
ulations of a large number of particle trajectories in the same flow 
discretized on a mesh. Numerical values of the parameters in these 
simulations: p v j pf = 1000, Fr = 0.1, t = 0.05. 


A playing the role of B and vice versa), and we readily ob¬ 
tain Aflg(t n ) = 2St/Fr for all n. This means that, for such 
a horizontal perturbation, only gravity matters and drives the 
particles towards the lower half-cell during their motion from 
B to A. 

Because we want to determine the critical Stokes number 
of inertial particles advected by numerical velocity fields, we 
have performed numerical computations of particle trajecto¬ 
ries using a discretization of the analytical velocity field (14). 
Typical clouds with various diameters are shown in Fig. 4. 
We then have checked whether or not particles have an ir¬ 
regular dynamics around the dividing streamlines, for various 
w’s between 0 and 12, and recorded the corresponding critical 
Stokes number below which this happened. This St c has then 
been compared to the theoretical expression (19) in Fig. 3. 
We observe that the agreement is good, so that the asymptotic 
analysis can be applied to numerical velocity fields also. 

The application of the asymptotic theory to numerical 
fields requires that the positions of points A and B, and the 
curvature-weighted integral (|uo| 2 ) c , be determined from the 
discrete velocity field. In the case of the cellular flow inves¬ 
tigated in this section, separatrix AB has been determined 
approximately by calculating a fluid point trajectory in the 
steady velocity field (e = 0), injected very close to point 
A. The curvature-weighted integral (|u 0 | 2 ) c is then com¬ 
puted. For our finest mesh we obtain (|uo| 2 ) c — —6.28, 
which agrees reasonably with the expected value —27r. 

Note however that the computation of (|u 0 | 2 ) c done in this 
section could also be performed from any velocity field ob¬ 
tained by solving the Navier-Stokes equations with a numer¬ 
ical solver. Then, the behaviour of particles near the separa¬ 
trix can be predicted by simply applying formula (11), or its 



FIG. 4: Test-case : typical clouds of particles transported in the 
vicinity of the unsteady cell flow (14), and initially injected near 
point A, with y > 0. Six particle diameters are used, corresponding 
to six different values of the Stokes number. All Stokes numbers are 
super critical for the lower separatrix (that is: Stf < St). Hence, 
particles who reach the lower circular separatrix AB are too iner¬ 
tial to re-enter the cell. Numerical values of the parameters in these 
simulations: p P /pf = 1000, Fr = 0.1, e = 0.05. 


dimensional counterpart (12), without doing any inertial par¬ 
ticle tracking. This is done in the next section, in the case of a 
cavity flow and of a backward-facing step. 


IV. APPLICATION TO NUMERICAL VELOCITY FIELDS 

To illustrate further these results and show their applica¬ 
bility to numerical solutions of the Navier-Stokes equations, 
we consider the classical examples sketched in Figs. 1 (Hi) 
and i{iv). We claim that the behaviour of particles can be 
analysed by calculating the various terms of Eq. (10), re¬ 
written here in dimensional form (keeping the same symbols 
for streamfunctions): 

— AV’o ~ 9 T p Xab — t p (|U 0 | 2 ) C - (20) 

Numerical method. For both the cavity and the backward 
facing step, numerical simulations using a finite-volume al¬ 
gorithm have been performed on a structured mesh using 
Gauss integration and linear interpolations. The prediction- 
correction PISO algorithm has been used 19 . Time stepping is 
done by means of a Crank-Nicolson scheme, and a CFL num¬ 
ber (Courant-Friedrichs-Lewy) below 0.1 is used in all cases. 
Particles are calculated by means of linear interpolation cou¬ 
pled to a Euler time-stepping method. They are allowed to 
slip on the walls, so that the separatrix streamline and the wall 
form a heteroclinic loop. For each geometry, two runs have 
been performed. First, a steady boundary condition has been 
used, together with a moderate Reynolds number, to obtain 
a steady velocity field with streamlines of the form of those 
sketched in Fig. 1. Then, the separatrix and the curvature- 
weighted integral (|XJo| 2 )c have been computed, and injected 
into formula (12). This equation also requires to determine 
the flux 6. This quantity, which is only flow dependent, can 
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FIG. 5: Critical diameter obtained from the asymptotic theory (solid 
line : cavity, dashed line : backward-facing step). Symbols indicate 
the diameters obtained by tracking a large number of inertial parti¬ 
cles (circles : cavity, triangles : backward-facing step). Numerical 
values of the parameters: ui = 26.2 rad/s, t = 0.05. 


be estimated after a run involving fluid points only, since the 
variation of streamfunction is also the flow rate between the 
two corresponding streamlines. In the vicinity of any point C 
of the separatrix we have: 

6 = max|AV’o(in)| ~ - tp 0 (N)\ 

where M and N denote the intersection between the border 
of the stochastic zone and the line normal to separatrix AB at 
C (see Fig. 2 (f)). A first order expansion of 'ipo(M) — ■ipo(N) 
leads to 

S ~ d|U 0 (C7)|, (21) 

where d = \AfN\ is the thickness of the stochastic layer at 
C. By measuring the thickness d from fluid point simula¬ 
tions, we estimate the flux of fluid points <5 from formula (21), 
and the critical response time r pc , or diameter I), from Eqs. 
(12) and (13). We have checked that the exact position of C 
on the separatrix does not affect the resulting S, as expected. 
These predictions have then been compared to those obtained 
by tracking a large number of particles injected near A , and 
by checking whether separatrix crossing occurred. 

Cavity flow. We first consider the flow above a two- 
dimensional open cavity sketched in Fig. 1 (in). This flow 
has been widely investigated, in various contexts 21 23 . The 
upstream flow has an imposed parabolic velocity profile with 
a cross-sectional average velocity = 0.13 m/s. The depth 
and width of the cavity are L = 0.01 m and W = 0.02 m re¬ 
spectively. The corresponding Reynolds number, based on the 
viscosity of air at 20°C, is Re = U/fL/iv ss 90, which corre¬ 
sponds to a moderate value. The flow is steady and character¬ 
ized by a separatrix streamline AB separating open stream¬ 


lines (above the cavity) from closed streamlines (in the cav¬ 
ity). We have chosen to introduce some unsteadiness by im¬ 
posing a pulsatile velocity upstream of the cavity: 

E4o(i) = £/^(l + esinwf) (22) 

with lu = 26.2 rad/s, which corresponds to a period of the or¬ 
der of the convective time over the cavity, and e = 0.05. The 
question of interest here is to determine under which condi¬ 
tions aerosols released outside (near point A) can be captured 
within the cavity, and remain there either permanently or tem¬ 
porarily. By computing a lagrangian fluid point trajectory near 
AB we get (|U 0 | 2 ) C ~ 2.5 10 -4 m 2 /s 2 < gX AB = gW ~ 
0.2 m 2 /s 2 , so that gravity dominates and Eq. (20) reads 

- AV’o -gr p X AB . (23) 

We therefore conclude that, in the steady case (AipQ (t n ) = 0 
for all n), the structure of the manifolds is similar to the one 
sketched in Fig. 2(a), and particles injected sufficiently close 
to A drop regularly in the cell and cannot exit. This behaviour 
persists in the unsteady case provided T P > T pc given by 


In contrast, particles such that T p < Tp C , injected close enough 
to the separatrix, will have a chaotic dynamics and wander in 
and out. The flux of fluid points along the stochastic zone <5 
has been obtained by measuring the thickness d of the stochas¬ 
tic zone around the separatrix AB, and using formula (21). 
We obtain 6 ~ 5.010” 5 m 2 /s. This leads to a critical re¬ 
sponse time r pc ~ 2.5 10” 4 s. To check these predictions we 
have computed inertial particle trajectories with varying den¬ 
sity ratios p p /pf. We then have controlled whether particles 
can exit the cavity or not. The critical diameter below which 
this happens has been plotted versus the density ratio (Fig. 
5 (circles)). The analytical D c (obtained from Eqs. (24) and 
(13)) has been plotted on the same graph (solid line), we ob¬ 
serve that the agreement is satisfactory. The systematic over¬ 
estimation of the theoretical diameter, compared to the numer¬ 
ical one, might be due to the numerical determination of the 
thickness of the stochastic layer d. Indeed, by measuring d 
from fluid points simulations, we slightly overestimate the ac¬ 
tual thickness. Therefore, the theoretical critical diameter is 
overestimated too. Snapshots of particle clouds with diame¬ 
ters D € [1, 50] microns are shown in Fig. 6 in the case where 
p p = 100 kg/m 3 . All particles are injected near A. The most 
inertial particles (in red), for which D > D c = 26.5 microns, 
drop in the cavity and remain there permanently. 

Backward-facing step. A similar analysis has been 
done with the backward-facing step flow sketched in Fig. 
1 (iv). This flow has been widely investigated numerically or 
experimentally 24-27 . It is an interesting configuration where 
centrifugal and gravity effects act in opposite directions if the 
step is downward. Indeed, we have gX AB > 0 and, be¬ 
cause the dividing streamline is clockwise, (|Uo | 2 )c < o. In 
our simulations, the upstream flow is horizontal with a cross- 
sectional average velocity U/f = 0.17 m/s. The vertical size 
of the step is L = 0.01 m and the Reynolds number is Re = 
U/^L/v w 115. Here also, for such a moderate Reynolds 
number, the flow is steady. Also, X AB = 0.057 m. The grav- 
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FIG. 6: Typical clouds of particles moving in the vicinity of the 
cavity flow. The thick line is the dividing streamline AB of the steady 
flow. Colours indicate particles’ diameter D. in m. All particles are 
injected near A, and their density is p p = 100 kg/m 3 . Those for 
which D > D c = 26.5 microns have a regular dynamics and cannot 
exit. Other parameters: U ^ = 0.13 m/s, cavity size 1 cm x 2 cm, 
lo = 26.2 rad/s, e = 0.05. 



FIG. 7: Typical clouds of particles moving in the vicinity of the 
backward facing step. The thick line is the dividing streamline AB 
of the steady flow. Colours indicate particles’ diameter D, in m. All 
particles are injected near A, and their density is p p = 100 kg/m 3 . 
Those for which D > D c = 34.5 microns (red particles) drop 
in the cell and are trapped permanently. They have a regular dy¬ 
namics. Those for which D < D c (blue) have a chaotic dynam¬ 
ics near the heteroclinic cycle and can exit in a complex manner. 
Other parameters: U^, = 0.17m/s, vertical step size L = 0.01m, 
cj = 25.2 rad/s, e = 0.05. 

ity term appearing in Eq. (20) is gXAB — 0.56 m 1 2 3 4 /s 2 and 
the curvature term is (|Uo| 2 ) c — —0.00145 m 2 /s 2 . Gravity 
therefore dominates the dynamics of particles near the sep- 
aratrix: in the steady case the structure of the manifolds is 
similar to the one sketched in Fig. 2(c), and particles injected 
very close to A drop regularly in the cell and cannot exit. This 
behaviour persists in the unsteady case provided T p > Tpc 
given by Eq. (12). In these simulations also unsteadiness is 
forced by imposing a pulsatile inlet velocity (Eq. (22)). The 
corresponding critical diameter is plotted in Fig. 5 (triangles). 
The flux 6 has been obtained by computing fluid points tra¬ 
jectories near the separatrix, and measuring the thickness of 


the stochastic zone. Then, by using formula (21) we get 
5 ~ 2.5 1 () -' m 2 /s. The corresponding critical diameter is 
shown on the same graph (dashed line). Fig. 7 shows particle 
clouds with various diameters between 32 and 38 microns. 

V. CONCLUSION 

The analysis presented in this study faithfully predicts the 
behaviour of low-Reynolds number heavy particles in the 
vicinity of the boundary of recirculation cells. It concerns 
the limit of vanishing Stokes numbers, that is particles with 
a finite but small inertia. For such particles, separatrix cross¬ 
ing depends on the competition between three mechanisms, 
namely flow unsteadiness, streamline curvature, and gravity. 
The former effect always contributes to separatrix crossing 
and chaotic dynamics (if the separatrix is part of a homoclinic 
or heteroclinic cycle). The other effects can either contribute 
to separatrix crossing or be opposed to it, according to the ge¬ 
ometry of the unperturbed flow. The criterion which has been 
derived in the present study enables to predict the various situ¬ 
ations. It has been first confronted to numerical results involv¬ 
ing an analytical flow discretized on a grid. Then, the criterion 
has been applied to numerical solutions of the Navier-Stokes 
equations, and we obtained satisfactory results. 

Note that the method can be applied to a wide variety of 
two-dimensional unsteady flows, even if the unsteadiness is 
more complex than the periodic time-dependence considered 
here. In particular, the method should work also for multi¬ 
modal perturbations, or even turbulent flows, provided the tur¬ 
bulence intensity is weak (i.e. the typical scale v! of velocity 
fluctuations is much smaller than the scale U of the average 
velocity). The only contribution of the flow unsteadiness to 
separatrix crossing is contained in the flux 5, which can be 
measured accurately from simulations of fluid point trajecto¬ 
ries and Poincare sections near the separatrix. 

The generalization of the asymptotic analysis to particles 
that are not much heavier than the fluid is the next step of this 
work. In this case additional contributions of the hydrody¬ 
namic force must be taken into account: pressure gradient of 
the undisturbed flow, added mass, history force, drag and lift 
due to the inertia of the displaced fluid. Even though the exact 
expression of the last two forces is still under debate, they are 
known to affect the complex dynamics of the particles (see for 
example Guseva el al/' for a detailed analysis of the effect of 
the history force). The generalization of the asymptotic analy¬ 
sis to such particles would be useful for the prediction of tiny 
objects in water flows. 
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